High-Throughput Sequencing of Oral Microbiota in Candida Carriage Sjögren’s Syndrome Patients: A Pilot Cross-Sectional Study

Background: This study sought to characterize the saliva microbiota of Candida carriage Sjögren’s syndrome (SS) patients compared to oral candidiasis and healthy patients by high-throughput sequencing. Methods: Fifteen patients were included, with five Candida carriage SS patients (decayed, missing, and filled teeth (DMFT) score 22), five oral candidiasis patients (DMFT score 17), and five caries active healthy patients (DMFT score 14). Bacterial 16S rRNA was extracted from rinsed whole saliva. PCR amplification generated DNA amplicons of the V3–V4 hypervariable region, which were sequenced on an Illumina HiSeq 2500 sequencing platform and compared and aligned to the SILVA database. Taxonomy abundance and community structure diversity was analyzed using Mothur software v1.40.0. Results: A total of 1016/1298/1085 operational taxonomic units (OTUs) were obtained from SS patients/oral candidiasis patient/healthy patients. Treponema, Lactobacillus, Streptococcus, Selenomonas, and Veillonella were the primary genera in the three groups. The most abundant significantly mutative taxonomy (OTU001) was Veillonella parvula. Microbial diversity (alpha diversity and beta diversity) was significantly increased in SS patients. ANOSIM analyses revealed significantly different microbial compositional heterogeneity in SS patients compared to oral candidiasis and healthy patients. Conclusion: Microbial dysbiosis differs significantly in SS patients independent of oral Candida carriage and DMFT.


Introduction
Sjögren's syndrome (SS) is a systemic autoimmune exocrinopathy characterized by mouth and eye dryness, fatigue, joint pain, and other symptoms [1]. As an effect of pathology, the oral manifestations of SS include xerostomia and hyposalivation, dental caries, oral candidiasis, dental erosion and/or abrasion, and others [2]. In 2015, untreated dental caries affected an estimated 2.6 billion adults globally. That year, oral conditions contributed to the declined health in 35 of 39 categories of cancer [3]. Among all adults with dental caries, patients with SS showed a significantly higher dental caries rate compared to non-SS patients [4]. SS patients are also more likely to have had teeth extractions [5] and have a higher rate of decayed, missing, and filled teeth (DMFT) [6].
The development of dental caries depends on the balance of pathological and protective factors. Pathological factors include acidogenic bacteria, inhibition of salivary function, and frequency of ingestion of carbohydrates [7]. Hyposalivation (defined as a whole salivary flow rate <0.1 mL/min without stimulation) is an important influence during the development of dental caries in SS [8]. Hyposalivation might delay the clearance of sugar

Participants
Because pSS has a strong female propensity, all participants in this study were female [29]. They were recruited at Peking University School and Hospital of Stomatology, Beijing, China, from January to August 2017 consecutively. The individuals in the Test group were diagnosed with primary SS according to the revised international classification criteria proposed by the American-European Consensus Group, 2002 [30]. They also met the 2016 ACR-EULAR Classification Criteria for primary Sjögren's Syndrome, with a total score ≥4 [31]. All participants were interviewed. Those who had not received antibiotics in the preceding 3 months and had no systemic diseases (including hepatitis, tuberculosis, and diabetes) and no smoking or alcohol abuse were enrolled. All enrolled participants provided signed informed consent. This study was approved by the Institutional Review Board of Peking University School and Hospital of Stomatology (#PKUSSIRB-201947088).

Clinical Examinations
A comprehensive oral examination was performed by an experienced dentist through visual/tactile methods after air drying of the mouth. The total numbers of teeth present, decayed (D), missed (M), filled (F) teeth (T)/surfaces (S), incisal/cusp caries, and root surfaces caries were recorded separately. The dental examination was completed within 15 min, and the findings were recorded by an assistant. The examiner checked the record after the examination. Individuals with more than 20 teeth were enrolled. To compare the oral health status, the inclusion criteria for SS patients (Test group) were DMFS ≥15, teeth with active caries ≥5, and unstimulated whole salivary flow rate (UWSF) <0.1 mL/min. The respective values for oral candidiasis patients (positive control group, Group P) and healthy control group (Group N) were ≥10, ≥1, and ≥0.1 mL/min. Patients with severe chronic periodontitis (probing depth ≥ 6 mm) were excluded. Concentrated oral rinse fluid was cultured on Sabouraud dextrose agar. Candida counts >200 colony-forming units (CFU)/mL indicated Candida carriage for SS and candidiasis patients. Finally, the total of 15 patients included 5 Candida carriage SS patients in the Test group (A01 to A05), 5 patients in Group P (B01 to B05), and 5 caries active healthy patients in Group N (C01 to C05). Details of the original data from clinical examination was uploaded as Supplementary Materials.

Whole Saliva Collection
All included patients stopped oral hygiene 24 h prior to sample collection and did not eat on the morning of the examination day. The participants were instructed to rinse their mouths thoroughly for 1 min with 10 mL of sterile saline. The saline was completely spit into a sterile graduated test tube. Each tube was stored at −80 • C until analyzed.

DNA Isolation, PCR Amplification, and High-Throughput Sequencing
All samples were washed twice with TE buffer (10 mM Tris-HCl, 1M EDTA, pH 8.0). Bacterial DNA was isolated and purified using a Wizard Genomic DNA Purification Kit (Promega, Madison, WI, USA) according to the manufacturer's instructions. The concentration of the purified DNA was determined using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) at an optical density of 260 nm (OD260). DNA integrity was verified by 1% agarose gel electrophoresis. For each individual sample, a standard concentration of 10 ng/uL DNA was diluted for the subsequent PCR assay and then stored at −80 • C.
High-throughput sequencing of all samples was performed by Cnkingbio Lab (Beijing, China). The V3−V4 hypervariable regions of the bacterial 16S rRNA gene were amplified by PCR using a GeneAmp 9700 device (Applied Biosystems, Foster City, CA, USA) from the isolated DNA, with primers 515F (GTGCCAGCMGCCGCGGTAA) and 806R (GGAC-TACHVGGGTWTCTAAT) and a unique barcode. PCR involved 10 min at 95 • C, followed by 35 cycles of 10 s at 95 • C, 30 s at 55 • C, and 30 s at 72 • C, and a final 10 min extension at 72 • C. The amplified PCR products were purified and checked with 2% agarose gel electrophoresis. The final amplicon preparation products were pooled in an emulsion PCR performed using a GS Lib-L kit (Roche Diagnostics GmbH, Mannheim, Germany) at a molecules-per-bead ratio of 0.7 according to the manufacturer's instructions. Sequencing was performed on a HiSeq 2500 system (Illumina, San Diego, CA, USA) using paired-end and dual-end sequencing for the V3−V4 region.

Data Processing and Statistical Analyses
The generated sequences from Illumina HiSeq 2500 system were annotated as unique target sequences (raw tags) by FLASH (V1.2.7, http://ccb.jhu.edu/software/FLASH accessed on 9 October 2017). Clean tags were generated by filter raw tags using QIIME quantity control process (V1.9.0, http://qiime.org/scripts/split_libraries_fastq.html, accessed on 9 October 2017) [32]. The quality filtering strategy was performed by tags truncation and tags length. The effective tags sequences were selected by redundancy removal to remove chimera and impurity sequences using Mothur software v1.35.1 (http://www.mothur.org/wiki/Download_mothur, accessed on 9 October 2017) [33]. The effective unique tag sequences pre-clustered using the single linkage pre-clustering approach with 97% similarity were defined as one operational taxonomic unit (OTU). Species annotation of OTUs were classified to study the species composition diversity of the samples. The representative sequences from each OTU were compared and aligned to the SILVA database (Version 1123, http://www.mothur.org/wiki/Silva_reference_files, Bremen, Germany, accessed on 9 October 2017).
Average age was compared by independent-samples t-test. Clinical data (DMFT, DMFS, incisal caries) were compared by the chi-square test in the three groups via SPSS 20.0 (IBM, Armonk, NY, USA). OTU differences across groups were analyzed by Metastat analysis using Mothur software v1.35.1 (http://www.mothur.org/wiki/Download_mothur, accessed on 9 October 2017). Alpha diversity was defined as the mean species diversity in sites (within-community) or habitats at a more local scale. Alpha diversity is the number of bacterial taxa and the proportional representation of each taxon per sample. The ACE, Chao, Shannon, and Simpson indices of alpha diversity were tested using the Wilcoxon rank-sum test. To examine the differences in the alpha diversity, data were analyzed using Bayesian methods by the limma package in R software (v2.2.0). Beta diversity was determined by principal component analysis (PCoA) based on weighted Unifrac matrix using R software (v2.2.0). Analysis of similarities (ANOSIM) determined the difference in bacterial composition and structure among the three groups.

Patients Characteristics
The baseline characteristics of SS patients and controls are presented in Table 1. With 15 participants, significant differences were not evident for the mean age, DMFT, and DMFS (p > 0.05) among the three groups (Table 1). Incisal caries was examined only in the Test group of SS patients.

General Outline
An average of 132,153 (ranging from 13,112 to 298,785) total pairs reads was obtained. An average of 129,221 combined reads as effective tags was obtained after splicing using Flash software. A total of 3399 OTUs, comprising 1016 in the Test group, 1298 in Group P, and 1083 in Group N, were obtained from 15 samples. The samples comprised 10 phyla, 22 classes, 43 orders, 67 families, and 121 genera. The 121 genera comprised 85 genera in the Test group, 101 in Group P, and 95 in Group N. The ten most predominant genera constituting the sequences detected in saliva are presented in Figure 1. The five most predominant genera in the Test group/Group P/Group N were Treponema (10.2%/9.9%/10.2%), Lactobacillus (9.5%/10%/8.4%), Streptococcus (7.6%/7.9%/8.5%), Selenomonas (8.6%/7.1%/7.8%), and Veillonella (7.9%/7.1%/7.8%). The numbers and proportions of sequences identified at the genus level were comparable among the three groups (p > 0.05).
Rarefaction curves of samples are presented in Figure 2. Significant differences (p < 0.05) in the abundance of the taxa from saliva was estimated at the species level using Metastat analysis. Compared with Group P, 19 OTUs were significantly different in the Test group, with 7 OTUs increased and 12 OTUs decreased. Compared with Group N, 26 OTUs were significantly different in the Test group, with only 1 OTU (OTU0001, Veillonella parvula) increased and 25 OTUs decreased. Low taxonomy abundance is defined as an abundance less than 1% [34]. More details about taxonomy abundance greater than 0.001 are presented in Figure 3 and Table 2. Rarefaction curves of samples are presented in Figure 2. Significant differences (p < 0.05) in the abundance of the taxa from saliva was estimated at the species level using Metastat analysis. Compared with Group P, 19 OTUs were significantly different in the Test group, with 7 OTUs increased and 12 OTUs decreased. Compared with Group N, 26 OTUs were significantly different in the Test group, with only 1 OTU (OTU0001, Veillonella parvula) increased and 25 OTUs decreased. Low taxonomy abundance is defined as an abundance less than 1% [34]. More details about taxonomy abundance greater than 0.001 are presented in Figure 3 and Table 2.

Biodiversity of the Salivary Microbiota
For alpha diversity, microbiota community richness was analyzed by the ACE and Chao indices, and the community diversity was analyzed by the Shannon and Simpson indices. For all four indices, significant differences were evident among the three groups (p < 0.05), as in Figure 4. Results of the beta diversity of the salivary microbiota determined by PCoA is shown in Figure 5. A clear separation of sample distribution was evident in the Test group, with roughly greater values than those from both Group P and Group N (principal components of 12.80% and 14.80%, respectively).

Biodiversity of the Salivary Microbiota
For alpha diversity, microbiota community richness was analyzed by the ACE and Chao indices, and the community diversity was analyzed by the Shannon and Simpson indices. For all four indices, significant differences were evident among the three groups (p < 0.05), as in Figure 4. Results of the beta diversity of the salivary microbiota determined by PCoA is shown in Figure 5. A clear separation of sample distribution was evident in the Test group, with roughly greater values than those from both Group P and Group N (principal components of 12.80% and 14.80%, respectively).   Similarity in bacterial composition is indicated by the distance between dots, with smaller distance denoting increased similarity. The spatial distribution of plots from Group P and N subjects overlapped, suggesting a similar set of genera. Plots from the Test group of subjects formed a cluster, having roughly greater values than the plot clusters from Group P and N subjects. Test group samples displayed a trend of dispersion along the PCoA2 axis. Group P and N samples displayed a trend of dispersion along the PCoA1 axis.
ANOSIM analysis also revealed significant differences in bacterial composition and structure among three groups (R = 0.304889, p = 0.005 < 0.01, respectively). Significant differences were evident in bacterial composition and structure between the Test group and Group P (R = 0.308, P = 0.043 < 0.05, respectively) and the Test group and Group N (R = 0.612, P = 0.003 < 0.01, respectively). No significant differences in bacterial composition and structure were evident between Group P and Group N (R = −0.032, p = 0.551 > 0.05, respectively).   Similarity in bacterial composition is indicated by the distance between dots, with smaller distance denoting increased similarity. The spatial distribution of plots from Group P and N subjects overlapped, suggesting a similar set of genera. Plots from the Test group of subjects formed a cluster, having roughly greater values than the plot clusters from Group P and N subjects. Test group samples displayed a trend of dispersion along the PCoA2 axis. Group P and N samples displayed a trend of dispersion along the PCoA1 axis.
ANOSIM analysis also revealed significant differences in bacterial composition and structure among three groups (R = 0.304889, p = 0.005 < 0.01, respectively). Significant differences were evident in bacterial composition and structure between the Test group and Group P (R = 0.308, P = 0.043 < 0.05, respectively) and the Test group and Group N (R = 0.612, P = 0.003 < 0.01, respectively). No significant differences in bacterial composition and structure were evident between Group P and Group N (R = −0.032, p = 0.551 > 0.05, respectively).

Figure 5.
Beta diversity evident in principal coordinate analysis (PCoA) based on weighted Unifrac matrix. Similarity in bacterial composition is indicated by the distance between dots, with smaller distance denoting increased similarity. The spatial distribution of plots from Group P and N subjects overlapped, suggesting a similar set of genera. Plots from the Test group of subjects formed a cluster, having roughly greater values than the plot clusters from Group P and N subjects. Test group samples displayed a trend of dispersion along the PCoA2 axis. Group P and N samples displayed a trend of dispersion along the PCoA1 axis. ANOSIM analysis also revealed significant differences in bacterial composition and structure among three groups (R = 0.304889, p = 0.005 < 0.01, respectively). Significant differences were evident in bacterial composition and structure between the Test group and Group P (R = 0.308, P = 0.043 < 0.05, respectively) and the Test group and Group N (R = 0.612, P = 0.003 < 0.01, respectively). No significant differences in bacterial composition and structure were evident between Group P and Group N (R = −0.032, p = 0.551 > 0.05, respectively).

Discussion
Our study findings strongly support the existence of distinct differences in the composition of the salivary microbiota in SS patients, independent of Candida carriage and DMFT. These findings further suggest that dysbiosis in oral microbiota might be a characteristic in SS patients.
Treponema, Lactobacillus, Streptococcus, Selenomonas, and Veillonella were the dominant bacterial genera in saliva. Treponema, especially T. denticola, is commonly associated with periodontal disease [35]. Lactobacillus, Streptococcus, Selenomonas, and Veillonella have been significantly associated with the progression of dental caries [20]. Lactobacillus and Streptococcus are well-known cariogenic bacteria. Yet, prior studies have reported discordant results. One study found a similar salivary abundance of Lactobacillus and Streptococcus compared to a healthy population [6]. In contrast, a higher abundance in SS patients than a healthy population was found for Lactobacillus in two studies [17], and a similar abundance for Streptococcus also in two studies [14] but low in another study [36]. The reported abundance of Selenomonas increased in three studies [6, 19,37] but was low in another [36]. Another predominant bacterial genera, Prevotella and Veillonella, was enriched in salivary microbiota in two studies [6,14]. The abundance of Neisseria was decreased in two studies [16,38].
In this current study, we did not observed differences in the abundance of any of these taxa. There are two possible explanations for this discrepancy. One is dental caries severity caused by oral cariogenic bacteria, and the other is the interaction between oral Candida and oral bacteria. Dental caries is a bacterial disease caused by dysbiosis in the complex microbiota community [21]. SS patients exhibit higher DMFT and DMFS. Our results agreed with the prior observations that with comparable DMFT, the salivary microbiota with most abundant genera of Prevotella, Veillonella, and Streptococcus in SS patients did not differ significantly from healthy controls [18]. Moreover, the relative abundance of Veillonella in SS patients with higher DMFT was reported to be higher than the abundance in healthy controls [6]. In this current study, the nonsignificant differences in the abundance of these cariogenic bacterial genera among the groups could be attributed to the comparable DMFT, leading to similar microbiota community dysbiosis.
Candida is a major cause of the opportunistic disease. The higher prevalence of oral Candida has been associated with hyposalivation and adverse effects on oral health [39]. Increasing evidence is solidifying the interaction between oral Candida and oral cariogenic bacteria, especially C. albicans [40]. In early in vitro oral biofilms, the presence of C. albicans alters the bacterial microbiota, leading to the presence of strictly anaerobic bacteria when oxygen is abundant [26]. The presence of C. albicans would increase the abundance of Streptococcus (particularly S. mutans), certain Lactobacillus species, and salivary Veillonella and Prevotella [26,41] and also increase the cariogenic virulence of biofilms [42] and even serve as potential "keystone" components of oral biofilms [43]. Oral Candida are also correlated with the number of untreated DMFT [44]. It was reported that Candida-infected individuals had 2.3 times more dental caries than those without Candida in their oral cavity [45]. C. albicans promotes caries activity that is dependent on PHR2 in a mixed microbial consortium [46]. In this current study, with similar oral Candida carriage and comparable DMFT, the abundance of the salivary microbiota at the genus level was not significantly different between SS patients with controls.
At the species level, the significantly increasing abundance of V. parvula (OTU001) in SS patients than controls was observed in this present study. The presence of Veillonella in the salivary microbiota has been related to the occurrence of dental caries [47]. Veillonella is a gram-negative anaerobic coccus capable of fermenting organic acids such as malate and lactate. This may result in increasing abundance of the salivary microbiota in SS patients [48]. The survival of V. parvula can alter the physiology of S. mutans by changing the expression of genes coding for many proteins [49]. This could contribute to more prevalent carious lesions in SS subjects. In this current study, even with comparable DMFT, significantly more V. parvula was detected in SS patients. These findings demonstrate the connection between V. parvula and the autoimmune disease SS. Elevated dysbiosis of V. parvula can initiate the progression of SS [50]. Hence, V. parvula, as an immunomodulatory commensal bacterium, may serve as a unique microbial biomarker for SS. This suggestion is supported by the results in other studies [6,14,51]. Uniquely, this present study is the first confirmation of the dominant role of V. parvula. Detection of the abundance of V. parvula might benefit the diagnosis, treatment, and monitoring of SS.
Measurement of microbial community diversity can be critical for understanding the composition of the microbiota [52]. Different from previous studies [16,17,19,38], alpha diversity analyses of each sample revealed significantly higher diversity and more equal distribution of abundance in samples from SS patients than from control patients in this present study. Our results are supported by prior studies [6, 14,53], which also reported higher alpha diversity in SS patients compared to healthy controls. Concerning beta diversity, the differentiation among sites or habitats contrasts with previous study findings [54]. In this present study, we observed more diverse bacterial composition of the saliva microbiota from SS patients. Similar results were previously reported in SS patients compared to a control group [36]. Normally, the difference in the oral microbial diversity of patients might be related to xerogenic medications, oral health status, saliva flow rates, salivary pH, severity of dental caries, and oral Candida carriage. Compared with oral Candida carriage and dental caries severity, the significant differences in microbial diversity might be a result of the effects caused by SS on the oral environment. ANOSIM analysis indicated SS itself might be associated with oral microbial dysbiosis. This was also a conclusion previously [14].
Our collective results suggest that SS patients with Candida carriage have different diversity and abundance of microbiota composition compared to patients with only Candida carriage and healthy patients, irrespective of DMFT. It remains unclear whether these changes observed in microbiota are the cause or consequence of SS. Interaction between the oral microbiota and autoimmune diseases is complex. For SS, activated CD4+ T cells and B cells infiltrate into the salivary glands [55]. Salivary gland destruction changes in saliva protein composition, which leads to microbiota dysbiosis in the oral cavity [50]. Changes in the oral microbiota may also be a factor for the onset of systemic autoimmune diseases, including SS [50]. Autoimmune diseases might be influenced by oral microbiota through the activation of Toll-like receptors, molecular mimicry, epitope spread, and antigen persistence [56]. The molecular mimicry theory posits that immune cells respond to self-antigens sharing epitopes with foreign antigens from the microbiota. Animal studies have demonstrated that peptides from oral bacteria could activate Ro60-reactive T cells [57], followed by the activation of B cells to plasma cells that produce SS antigen A [56]. Future studies should examine interactions between microbiota imbalances and autoimmune responses in the development and progression of SS, especially the role of V. parvula.
The diagnosis and management of SS represent a challenge for the clinician, especially the dentist. Accompanying various lesions and pseudoconditions of the oral mucosa, the differential diagnosis must be made [58]. This study found significant dysbiosis in the oral microbiota of SS patients through the marked change in salivary microbiota compared with controls. This study provided a potential way to help with diagnosis and management. There are some limitations of this study. First, the sample size is small, with only five in each group. This may have limited the power to detect significant associations. However, for the high-throughput sequencing based experiments, power and sample size calculations are not established. In some oral microbiome-related references, comparisons have involved 11 subjects divided into 2 groups [59] and 10 SS patients divided into 2 subgroups [60]. These two studies attained reasonable results and conclusions. Therefore, the five subjects for each group in this study were likely sufficient to prove the composition diversity and identify the most important bacteria in samples [61]. Second, 16S rRNA gene sequencing remains the gold standard of environmental microbiology. However, species level classification is not reliable with 200-400 bp reads of the 16S rRNA gene, and analyses have been confined to genus and higher taxonomies [62]. Finally, as one of the clinical characteristics, the UWSF in primary SS patients was −0.43 mL/min compared to controls [63]. It is difficult to conclude whether the microbial dysbiosis was caused by SS or hyposalivation. Future studies should involve many more patients (females and males) and controls to provide a more comprehensive view of microbiota composition. Quantitative real-time PCR should be used with standard strains to determine the absolute abundance of each taxon.

Conclusions
Oral saliva microbiota changes are documented for the first time in oral Candida carriage SS patients compared to control patients of comparable age and DMFT. The present findings extended the depth of previous findings. V. parvula might be a unique immunomodulatory commensal bacterium and perhaps a microbial biomarker for SS. The abundance of dominant bacteria did not differ significantly between the groups. However, the distinct community structure diversity identified microbial dysbiosis as a significant characteristic in SS patients, irrespective of oral Candida carriage and dental caries status.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.